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A  Generalized  Finite  Element  Method  for 
Solving  the  Helmholtz  Equation  in  Two 
Dimensions  with  Minimal  Pollution 

Ivo  M.  Babuska  Prank  Ihlenburg  Ellen  T.  Paik 

Stefan  A.  Sauter 


Abstract 

When  using  the  Galerkin  FEM  for  solving  the  Helmholtz  equation  in 
two  dimensions,  the  error  of  the  corresponding  solution  differs  substantially 
from  the  error  of  the  best  approximation,  and  this  effect  increases  with 
higher  wave  number  k. 

In  this  paper  we  will  design  a  Generalized  Finite  Element  Method 
(GFEM)  for  the  Helmholtz  equation  such  that  the  pollution  effect  is  mini¬ 
mal. 

1.  Introduction 

Boundary  value  problems  governed  by  the  Helmholtz  equation  arise  in  many  phys¬ 
ical  applications,  as  for  example  the  scattering  of  a  wave  from  an  elastic  body. 

For  this  kind  of  problems  the  computational  domain  consists  typically  of  the  finite 
domain  of  the  elastic  body  coupled  with  the  unbounded  exterior  domain  for  the 

scattering  field.  In  the  unbounded  domain  the  scattered  wave  is  described  by  the  _ 

classical  Helmholtz  equation  with  the  Laplace  operator  as  the  principal  operator,  _ 

In  the  elastic  body  the  equation  is  of  the  same  so-called  Helmholtz  type  with  the 
Laplace  operator  replaced  by  the  elasticity  operator.  □ 

In  order  to  solve  numerically  such  a  coupled  scattering  problem,  the  finite -d  □ 

element  method  is  typically  applied  in  the  elastic  body.  For  the  approximation' - - - 

of  the  scattered  wave,  various  approaches  are  used,  such  as  the  boundary  element 
method,  the  method  of  infinite  elements  coupled  with  finite  elements  and  the  finite- . . 
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element  method  where  the  unbounded  exterior  domain  is  replaced  by  a  bounded 
artificial  domain  with  suitable  boundary  conditions  on  the  artificial  boundary  (see 
[8]). 

In  this  paper  we  will  address  the  Helmholtz  problem  for  the  Laplace  operator 
on  a  bounded  domain  as  the  model  problem  which  characterizes  the  behavior  of 
the  finite  element  method  for  both  the  Helmholtz  problem  of  elasticity  and  wave 
scattering. 

It  is  known  and  understood  (see  [5])  that  the  accuracy  of  the  Galerkin-FEM 
deteriorates  with  increasing  wave  number.  To  be  more  concrete  we  have  to  intro¬ 
duce  a  norm  to  measure  the  accuracy  of  our  FE  solution.  Let  e  be  the  accuracy 
we  would  like  to  achieve.  Then  there  exists  a  number  of  elements  no  =  no  (e) 
such  that,  in  the  corresponding  finite  element  space,  there  is  a  function  called 
“best  approximation”  with  an  error  less  than  or  equal  to  e.  Usually  the  FE  so¬ 
lution  needs  more  elements  to  get  the  same  accuracy  (say  rigai  (e)).  For  standard 
elliptic  problems,  the  Galerkin  method  is  quasi-optimal,  meaning  that  the  ratio 
Ugai/no  is  a  constant.  For  the  Helmholtz  equation  the  situation  is  different.  In 
this  case  the  ratio  Ugai/no  goes  to  infinity  with  increasing  wave  number.  We  call 
this  non-robust  behavior  with  respect  to  the  wave  number  the  “pollution  effect” . 

A  generalization  of  the  FEM  was  introduced  in  [1]  and  is  called  Generalized 
FEM  (GFEM).  This  method  covers  practically  all  modifications  of  the  FEM  which 
lead  to  a  sparse  system  matrbc.  In  [2]  two  of  the  authors  have  defined  a  GFEM 
called  stabilized  FEM  for  the  Helmholtz  equation  in  ID  with  the  property  that 
UstabUized/no  is  a  constant  independent  of  the  wave  number.  However,  in  the  same 
paper  it  was  proved  that  in  the  two-dimensional  case  there  exists  no  GFEM  such 
that  the  ratio  noFEM/no  is  bounded  with  respect  to  the  wave  number. 

To  explain  the  goal  of  this  paper  we  consider  the  discretization  of  the  Helmholtz 
equation  separately  from  the  discretization  and  incorporation  of  the  boundary 
conditions,  as  with  the  finite  difference  method.  In  our  paper  we  focus  on  the  ap¬ 
proximation  of  the  Helmholtz  equation  in  the  interior  of  the  domain  by  a  GFEM. 
The  approximation  of  the  DtN  mapping  for  the  definition  of  the  boundary  con¬ 
dition  and  the  effect  of  its  discretization  is  not  the  subject  of  our  investigation. 
In  matrix-algebraic  terms,  our  task  is  to  define  the  interior  stencil  of  the  system 
matrix  in  such  a  way  that  under  an  optimal  modeling  of  the  boundary  conditions 
the  ratio  ristabiiized/no  increases  as  slowly  as  possible. 

The  importance  of  a  proper  approximation  of  the  Helmholtz  operator  was 
worked  out  in  [2],  where  it  was  shown  that,  if  the  interior  stencils  lead  to  a 
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pollution,  this  effect  cannot  be  countered  by  any  discretization  of  the  boundary 
conditions. 

Our  paper  is  organized  as  follows: 

In  the  next  section  we  will  study  a  one-dimensional  model  problem.  We  will 
define  the  GFEM  for  the  Helmholtz  problem.  We  will  show  that  the  corresponding 
discretization  error  is  directly  related  to  the  difference  of  the  so-called  discrete 
wave  number  k  with  the  exact  one.  The  difference  k—k  is  called  “phase  lag  .  Using 
these  results,  we  are  able  to  construct  a  FEM  with  the  property  that  k  —  k  =  0 
which  additionally  satisfies  the  usual  consistency  conditions.  We  will  prove  that 
this  GFEM,  called  “stabilized  finite  element  method”  (SFEM)  has  no  pollution. 

In  Section  3  we  will  define  the  GFEM  in  two  dimensions  and  a  2-D  analogy 
to  the  phase  difference  k  -  k.  It  will  turn  out  that  in  2-D  every  GFEM  has  the 
phase  lag.  We  will  prove  that  this  phase  lag  leads  to  a  pollution  term  in  the  error 
estimates. 

In  Section  4  we  will  define  and  explain  a  measure  of  the  approximation  quality 
of  the  GFEM  discretization  for  the  Helmholtz  equation. 

In  Section  5  we  will  define  a  GFEM  called  Quasi-Stabilized  FEM  (QSFEM) 
which  leads  to  the  smallest  possible  pollution. 

In  the  last  section  we  will  present  the  results  of  a  2-D  implementation,  where 
we  compare  the  quality  of  the  quasi-stabilized  FEM  with  the  usual  Galerkin  FEM 
and  a  further  GFEM  called  generalized  least  squares  finite  element  method  (GLS- 
FEM).  The  latter  method  was  developed  by  Thompson  and  Pinsky  (rf.  [9])  based 
on  a  paper  of  Harari  and  Hughes  ([4]). 

2.  One-dimensional  model  problem 

In  this  section  we  will  consider  a  one-dimensional  model  problem  and  explain  why 
the  accuracy  of  the  Galerkin  FEM  deteriorates  with  increasing  wave  number  k. 
This  non-robust  behavior  with  respect  to  k  is  called  the  pollution  effect  and  will 
be  defined  formally  in  Definition  2.3.  We  will  explain  how  the  pollution  effect 
is  related  to  the  underlying  discretization  method.  By  using  that  investigation 
we  are  able  to  construct  a  so-called  generalized  finite  element  method  (GFEM) 
having  no  pollution. 

To  fix  the  notation  let  us  consider  the  following  one-dimensional  model  prob¬ 
lem 

-u''-eu  =  fmQ:={0,l)  (2.1) 
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with  boundary  conditions 


«(0)  =  0, 

iku  (1)  +  u'  (1)  =  0. 

To  apply  the  finite  element  method  to  this  equation  we  write  (2.1)  in  a  variational 
formulation,  seeking  u  G  V  :=  {u  G  (O)  |  v  (0)  =  0}  such  that 

a  (u,  v)  :=  f  u'v'  -  k^uvdx  +  iku  (1)  v  (1)  =  /  (w)  (2.2) 

Jo 

is  fulfilled  for  all  u  G  V.  Here  and  in  the  following  we  assume  that  /  lies  in  the 
dual  space  V'  =  (f^). 

Further,  let  {xi}o<i<„  denote  a  set  of  grid  points  0  =  Xo  <  Xi  <  . . .  <  a:„  =  1. 
The  finite  element  grid  r  consists  of  the  intervals  {[a:m-i,a:m]}i<m<n- 
size  h  is  defined  by 

h  max  {xm  -  Xm-i)  ■ 

l<m<n 

We  consider  here  only  the  h- version  of  finite  elements  and  define  Sh  as  the  space 
of  continuous  functions  which  are  linear  on  each  interval.  A  study  of  the  p- version 
of  the  Galerkin-FEM  for  the  Helmholtz  equation  can  be  found  in  [7]. 

2.1.  The  GFEM  for  the  Helmholtz  equation  in  ID 

The  Generalized  Finite  Element  Method  (GFEM)  was  first  introduced  by  Babuska 
and  Osborn  (rf.  [1]).  The  idea  is  to  introduce  local  mappings  which  transform 
the  usual  finite  element  basis  functions  to  another  local  basis.  In  the  mentioned 
paper  these  local  mappings  were  designed  in  such  a  way  that  the  resulting  method, 
applied  to  a  differential  equation  with  highly  non-smooth  coefficients,  converges 
with  optimal  rate. 

For  our  purpose  we  define  the  GFEM  in  the  following  algebraic  way.  The 
GFEM  is  a  method  which  defines  a  tri-diagonal  matrix  A  G  C”^"  and  a  linear 
mapping  Q  :  ->  C".  The  solution  of 

Au  =  b 

with  b  :=  Q  (/)  is  then  identified  with  a  finite  element  function  by  the  basis 
representation  ^ 

Ufe.  (^)  “  5!^  ip) 

m=i 
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with  the  usual  local  nodal  basis  of  Sh-  The  function  u/e  serves  as 

an  approximation  of  the  exact  solution  of  f2.2).  How  the  Galerkin  FEM  can  be 
written  in  this  notation  can  be  found  in  the  following 

Example  2.1.  The  Galerkin  FEM  is  characterized  by  the  tri-diagonal  matrix 

Ajj  “  a  <^i) 

and  the  mapping  Q  defined  by 

hi  :=  mm  :=/(^i). 

2.2.  Error  analysis  for  the  Galerkin  FEM 

To  measure  the  accuracy  of  the  FEl-solution  we  have  to  introduce  suitable  norms. 
We  will  consider  the  and  "H^-seminorm  defined  by 

Iloilo  “  /  « (^)  “  (3^)  ^ 

J  o 

and 

Ihll.  :=  ll'‘'llo- 

From  the  approximation  theory  it  is  well  known  that  for  every  function  u  € 
(fi)  there  exists  it/,  €  5/,  such  that 

\\u-u,\\^<Ch^-^\\uX 

for  j  G  {0, 1}.  The  dependency  of  the  relative  error  on  h  and  k  is  discussed  in  the 
following 

Theorem  2.2.  Let  the  right-hand  side  /  of  (2.2)  be  in  C?  (fi).  We  assume  that 
the  exact  solution  u  of  (2.2)  is  oscillating  in  the  sense  that  forO  <  s  <t  <2 

<  Ck*-*  (2.3) 

is  satisfied.  For  j  G  {0, 1},  the  best  approximation  uipt  G  Sn  with  respect  to  the 
W -seminorm  is  defined  by 

uipt  ■■=  arg  inf  Iju  -  Uh\\j  ■ 
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The  {auction  u^pt  satisfies 


^  -  <pt\\j 


<  C  [hkf-^ . 


Proof.  Using  (2.3)  we  obtain 


4>t  = 


u  —  u, 


opt 


n 


Remark  1.  From  the  Theorem  2.2  we  see  that  the  accuracy  of  the  optimal  ap¬ 
proximation  depends  only  on  the  number  of  the  elements  in  one  wavelength  of 
the  solution,  i.  e.,  depends  only  on  the  value  k  *  h.  To  relate  this  value  to  the 
accuracy  of  the  solution  is  a  “rule  of  the  thumb”  in  engineering  computations. 

We  say  that  the  GFEM  has  the  pollution  effect  if  it  is  possible  that  Copt  is 
small  but  the  error  of  the  GFEM  solution  is  arbitrarily  large.  The  details  are  in 
the  following 

Definition  2.3  (pollution  effect).  For  j  e  {0, 1},  let  the  error  of  the  GFEM- 
solution  u fe  be  defined  by 

i  h-Ufe\\j 

IWIIi  ' 

If  this  error  can  be  estimated  by 

e},<Ci(khf-^-i-C2k‘’  (khy. 

with  s  >  0  and  in  addition,  there  exists  right-hand  sides  for  problem  (2.1)  such 
that  the  corresponding  finite  element  error  can  be  estimated  from  below  by 

4e>C3k^ikhy 

with  r>0,  then  we  say  that  the  GFEM  has  the  pollution  effect. 
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In  view  of  the  Theorem  above  it  is  obvious  that  the  error  of  the  best  approx¬ 
imation  is  small  if  kh  is  small,  while  the  condition  s,  r  >  0  in  the  definition  of 
the  pollution  effect  has  the  consequence  that  for  sufficiently  large  k  the  quantity 
k*  {khY  can  be  arbitrarily  large,  i.e.  the  “rule  of  the  thumb”  mentioned  in  the 
Remark  1  does  not  lead  to  an  accurate  solution  if  A:  is  large. 

The  following  Theorem  shows  that  the  Galerkin  FEM  for  our  model  problem 
has  the  pollution  effect. 

Theorem  2.4.  Let  the  exact  solution  of  (2.2)  be  oscillating  in  the  sense  of  (2.3). 
Let  us  assume  that  our  grid  is  uniform,  which  means 

h>  —  Xffi  —  Xx  Xx—\  Vm,  %  G  {1, 2, ... ,  n) 

Let  the  right-hand  side  of  (2.2)  /  G  (fl)  and  Ugai  be  the  Galerkin  solution. 
Then  for  hk  <1  the  error  estimate 

el,  <  C  (kh)  +  C,k  (khf 

holds.  The  error  in  the  C^-norm  can  be  estimated  by 

eU  ~  <  (a + C,k)  (khf . 

Iloilo 

Proof.  In  the  proof  of  [5,  Theorem  5]  it  was  shown  that 

<(CiA  +  C2(tt)')  ||«"|| 

holds.  In  conjunction  with  (2.3)  we  obtain  the  desired  estimate  of 
The  £^-estimate  was  proven  in  [6,  Theorem  4]. 

■ 

Numerical  computations  in  [5]  shows  that  the  error  estimates  are  optimal,  i.e. 
there  are  cases  where  and  are  bounded  from  below  by  the  same  pollution 
term  as  from  above.  For  a  theoretical  investigation  see  Theorem  2.6  together  with 
Lemma  2.5. 
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2.3.  Relation  of  the  finite  element  error  to  the  discrete  wave  number 

In  this  section  we  will  explain  how  the  pollution  effect  is  related  to  the  difference 
of  a  discrete  wave  number  and  the  exact  one.  Later,  we  will  study  this  effect 
in  two  dimensions  as  well.  We  include  here  the  one-dimensional  investigation 
because  the  main  ideas  are  more  visible  than  in  two  dimensions.  To  avoid  too 
many  technicalities  we  consider  the  following  model  example  with  Robin  boundary 
conditions  on  both  sides  . 

-u"-k\  =  0  inf2  =  (0,l)  (2.4) 

—u'  (0)  —  iku  (0)  =  —2ik, 
u'  (1)  —  iku  (1)  =  0. 

If  not  stated  otherwise,  we  assume  throughout  this  chapter  that  fl  is  partitioned 
into  intervals  having  constant  length  h.  It  is  easy  to  check  that  the  exact  solution 
of  (2.4)  is  given  by 

u(a;)  =  e**^ 

We  consider  a  GFEM  of  the  form 

Du  =  b  (2.5) 

with  the  (n  +  1)  X  (n  -t- 1)  matrix  D 

Di  N 
N  D2  N 
N  D2  N 

N  •••  ••• 

D2  N 
N  Di 

and  the  right-hand  side  vector 

b  =  {-2ik,  0, 0, . . . ,  0)^ . 

We  assume  that  the  elements  of  the  matrix  D  can  be  expanded  in  a  Taylor  series 
of  the  form 

D,  =  h-'(l-  ikh  +  a„  (kh)^  +  i  (kh)  Y,  0^  (khf"  +  O  ((»)"”+“)) 

\  n=l  n=l  / 
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(2.6) 


Lh  =  fc-‘f2  +  ^7„(fcA)^  +  0((W)^' 

\  n=l 

N  =  h-^{-l  +  f,K{khf"  +  0{(khf'^^ 

\  n=l 

oti  +  A  =  -\  7i  +  25i  =  -1. 

These  assumptions  are  very  natural  in  view  of  the  underlying  equations.  Some 
comments  for  the  first  three  conditions  are  given  later  in  Remark  3.  The  impact 
of  the  last  two  equations  will  become  clear  in  the  proof  of  Theorem  2.6. 

The  system  (2.5)  can  be  solved  explicitly.  For  this  purpose,  we  will  use  the 
concept  of  the  discrete  Fourier  transform  to  solve  finite  difference  equations.  Al¬ 
ternatively,  one  could  employ  the  theory  of  fundamental  systems  for  our  one¬ 
dimensional  model  problem.  We  prefer  the  first  method,  because  the  discrete 
Fourier  transform  can  be  extended  straightforwardly  to  the  higher-dimensional 
case.  In  contrast  to  this,  the  theory  of  fundamental  systems  is  applicable  only 
in  the  one-dimensional  case  because  the  number  of  homogenous  solutions  of  a 
second-order  differential  equation  in  2-D  without  boundary  conditions  is  infinite. 

We  start  by  computing  the  discrete  symbol  of  a  difference  scheme.  For  this 
purpose  we  introduce  the  discrete  Fourier  transform  of  a  complex  vector  u  = 

{^m}„,ez={..., -1,0,1,...} 

a  (f)  :=  CFu)  (0  :=  f 

m=— oo 

For  a  difference  scheme  with  constant  coefficients  given  by 

p 

(Au)„  =  XI 

1--P 

the  discrete  Fourier  transform  can  be  computed  as 
(Ai)(f)  =  ESJ=-ooEL-pA,u„+,e‘^ 

=  Ai  A,e-“f  ES=_,o 

=  EL-p  A,e-«  ES=-oo  =  (Ef.-,  Aie-'-f)  u  (Q  . 
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The  function  a(^)  :=  (Ef=-p  is  called  the  discrete  symbol  of  the  difference 

operator  A.  Let  denote  the  zeros  of  a  in  the  interval  [-7r,7r[,  i.e. 

a(6)  =  0 

V/6  {-p, -p  +  l,...,p} 

e  [-7r,7r[. 

By  the  theory  of  finite  difference  schemes  it  follows  that  the  vector 

■■=  E  m  e  z 

l=-p 

satisfies 

At/  =  0. 

Simple  computations  yield  that  the  discrete  symbol  of  the  difference  operator 
which  corresponds  to  the  GFEM  (2.5)  is  given  by 

d(^)  =  D2  4-2Arcos^. 

The  zeros  of  the  symbol  are  given  by 

^  =  ifc 


i  =  iarcco.(-^).  (2.7) 

The  number  k  in  this  context  is  called  the  “discrete  wave  number” .  Consequently, 

satisfies  equation  (2.5)  for  indices  2  <  j  <n  —  1,  i.e. 

(Du)j.  =  0  Vy  €{2,3,...,n-l}. 

The  constants  Ci  and  C2  are  determined  by  the  equations 

(Du)i  =  —2ik, 


(Du)„  =  0, 


10 


in  other  words,  by  the  boundary  conditions.  They  are  given  explicitly  by 

“  0?ainfc+2DiArain(i(l-A))+N*8iii(fc(l-2/.)) 

(2.8) 

=  (p2  8mjfc+2PiArsin{jk(l-/»))+Af*8iii(*(l-2/»)) 

with  h  :=  1/n. 

To  summarize  the  explanations  above,  we  state  that  the  solution  of  our  tri¬ 
diagonal  system  of  linear  equations  (2.5)  is  given  by 

uj  =  Ci/^^  +  C2e-^^^  0<j<n,  (2.9) 

where  the  discrete  wave  number  k  is  given  by  the  zeros  of  the  discrete  symbol  of 
the  underlying  difference  operator  and  is  independent  of  the  boundary  conditions. 
The  boundary  conditions  then  determine  the  constants  Ci  and  C2- 

Now,  we  will  study  the  relative  error  of  the  GFE-solution  in  the  £^-norm  given 


The  functions  <j)j  are  the  usual  linear  nodal  basis  of  the  finite  element  space  Sh¬ 
it  turns  out  that  the  error  eo  is  directly  related  to  the  distance  of  the  wave 
number  k  from  the  discrete  wave  number  k.  Therefore,  before  we  start  to  estimate 
eo,  we  will  estimate  k-k.  The  details  are  given  in  the  following 

Lemma  2.5.  Let  kh  be  bounded  and  conditions  (2.6)  be  satisfied.  Then  either 

k  —  k  =  0 

or  there  exist  constants  qo,  qi  independent  of  k  and  h  but  possibly  dependent  on 
7  and  6  of  (2.6)  such  that 

qok  (khy^  <  |A:  -  ^1  <  qik  {khy° 

with  So  >  2. 

For  the  Galerkin-FEM  the  estimate  above  holds  with  Sq  =  2. 
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Proof.  The  discrete  wave  number  was  given  by  (2.7).  In  view  of  that  equation 
we  compute  the  quotient  Using  (2.6)  we  get 

2  +  ^  1  ^  /  I  ^ +  f  ,,  (tt)- . 

2N  2(-l+E“,4(*A)^)  V2'  r  ’  tSi 

The  modified  wave  number,  given  by 


i  arccos  +  Q71  + 


{khf  +  £  {khy 

$=2 


can  be  expanded  about  kh  =  0  ^s 


=  I  (\/{-7i-2i.)tA  +  E  (“)' 


Now  it  is  clear  why  we  imposed  the  fourth  condition  in  (2.6).  Under  this  assump¬ 
tion,  the  term  =  1  and  the  discrete  wave  number  converges  towards 

k  as  h-^Q.  The  equation  above  can  then  be  written  in  the  form 

k-k  =  0{k  {khy'>) 

with  even  Sq  ^  2.  Only  in  the  case  of  t*  =  0  for  all  s  >  1  we  obtain 

k  —  k  =  0. 

In  [5]  it  was  shown  that  for  the  Galerkin  method  k  =  k  +  O  ^k  {kh)  ^  holds  if  kh 
is  bounded.  ■ 

In  the  sequel  the  number  Sq  is  given  by  the  Lemma  above. 

Now,  we  will  estimate  the  relative  £^-error  eo  from  above  and  below.  The 
details  are  in  the  following 


Theorem  2.6.  Let  us  assume  that  the  considered  GFEM  has  the  property  k^k. 
This  means  that  Sq,  defined  above  is  Snite.  Let  kh  and  k  (fch)*®  be  bounded  and 
k  sufficiently  large. 

Then,  the  error  eo  of  the  GFE-approximation  of  the  solution  of  (2.4),  namely 
e**'®,  can  be  estimated  by 

c\k-k\<eo<C  {khf  +  c\k-k\. 
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Proof.  For  the  following  analysis,  we  will  use  the  interpolant  u'f^t  of  func¬ 
tion  e*"®,  i.e. 


Using  the  fact  that  ||e**'®||  =  1,  the  error  bq  can  be  estimated  from  above  by 
eo  =  II  e***  -  -H  e**®  -  -f  -  Ej=o  (a^)  | 


<  lie**  -  e'**!  +  Ije'*"  -  uL||  +  ||«L  -  D”=(,  Uj^j  (a!)|| 


(2.10) 


We  will  now  estimate  the  three  terms  on  the  right-hand  side  above  separately. 
Considering  the  first  term  we  get: 


||e»*  -  (e*'  -  e“*)  (e-““  -  e"^)  dx 


=  2(l-i 


(2.11) 


8m(Jb— i) 
k-k 


By  our  assumption  we  know  that  k  —  k  is  bounded,  thus  the  Taylor  expansion 
about  k  —  k  =  0  results  in 


jgifcx  _  gifcx||  <  -I-  O  (|fc  -  kf^ 


The  second  term  on  the  right-hand  side  of  (2.10)  can  be  estimated  by  using  a 
standard  interpolation  argument 

||«“'  -  «L||  <  Ch^  ||(e‘^)"||  =  C  {hif . 

Using  k  =  k  +  0{k  {khY°)  (see  Lemma  2.5),  we  obtain 

|e*-«Ll<C(«)^ 

For  the  last  term  in  (2.10)  we  proceed  as  follows: 


=  E  (4.t-u),M,,,n(uL-u), 
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with  the  mass  matrix  M  defined  by 


<j)i  (j;)  <j>j  (x)  dx 


and  (u^nt)  :=  (hm).  It  is  well-known  that  if  the  FE-space  corresponds  to  a 

uniform  grid,  the  jC^-norm  is  equivalent  to  the  weighted  Euclidean  norm,  resulting 

“  n  « 

“L  -  W  2  *  E  |(“L  -  ")„f  =•  ll"“  "  Al  ■ 

j=0  tn=0 

The  norm  |Hlp  is  called  the  /^-norm.  For  a  vector  u,  we  introduce  the  convention 

llUm|&  \Umf  . 

m=0 

Using  the  definition  of  and  (2.9)  we  obtain 

“L  -  W  ^  |k“"*  -  -  C2e-*i\  <  |1  -  C.|  ||e“-’'‘||,,+|C2|  . 

j=0 

Direct  calculations  yield 

e‘^^'  \  =  hf2  =  h  1  =  (n  -h  1)  h  <  2, 

*  j=o  j=o 


4nt  -  (^)  <  2  (|1  -  Cl  I  -t- 1(72 1) 

j=0 


Let  us  first  estimate  the  constant  C2  of  (2.8): 

-fce**  (Di  +  Are-**") 

a=  - ^ ^ 


Dl  sin  k  +  2Di  AT  sin  (k  (1  -  h))  -h  sin  {k  (1  -  2h))  ‘ 


Replacing  fc  by  A:  +  e  and  inserting  the  formulae  for  Di  and  N  (cf.  2.6)  into  the 
definition  of  C2  and  expanding  C2  as  a  Taylor  series  about  e  =  0  and  AA:  =  0,we 
get 

(72=  l  +  2(ai+A)  kh  +  Cikhf. 


kh  +  C{khf. 
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In  view  of  this  representation  it  is  clear  why  we  imposed  the  fourth  condition  in 
(2.6).  The  constant  then  can  be  estimated  by 

C2<C{khf.  (2.12) 

Applied  to  |1  -  Cl  I,  the  same  arguments  results  in 

|l-Ci|  <C(jfch)".  (2.13) 

Combining  all  estimates  above  we  conclude  that  the  inequality 

(Cie^o  +  Cje-**)  (i)  <  C  (khf  +  C\k-k 

i=o  II 

is  fulfilled,  completing  the  proof  of  the  upper  estimate. 

For  the  lower  estimate  we  proceed  in  the  following  way.  The  error  eo  can  be 
written  in  the  form 

eo  =  -  Cie‘**  -  026“***)  +  ^Cie‘*®  +  Cze'**®  -  ^  (x)  j  . 

By  definition  we  have  that  (^)  ^  interpolant  of  Cie**®  +  C2e  ***. 

Therefore,  we  know  that 

+  C2e-‘*®  -  ^  (x)  <  Ch^  ||(Cie**®  +  C2e-‘*®)1  <  C  {khf  (|Ci|  +  IC2I) . 

i=o  "  ' 

Using  (2.12,2.13)  we  obtain  that 

Cie‘*®  +  C2e-‘*®  -  S  (x)  <  {khf . 

i=o  11 

Thus,  eo  can  be  estimated  by 

eo  >  e**®  -  Cie‘*®  -  C2e“**®  -Cm  {khf  >  min  le**'®  -  ue**®  -  ve“**®  -Cm  {khf . 

Let  L  :=  •  y  where  [gj  denotes  the  largest  integer  less  than  or  equal  to  q. 

Let  II -11^  be  defined  by 

u{x)u{x)dx. 
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Using  this  norm,  eo  can  be  estimated  by 


>  min  e**®  -  we**®  -  ve  **®  ^  -  Cm  (kh) 


u,v6C 


The  computation  of  this  minimum  is  easy  in  view  of  the  orthogonality  of  e**®  and 
g-ifc®  (0,  L).  Introducing  e  =  k-k,  we  get  that  the  minimum  is  achieved  for 


=  1  =  ^ 

L  Jo  teL 

=  11’-  =  f"-]  . 

L  Jo  i(2k  +  e^  L  i  (2k  +  €^  L 


We  assumed  that  k  is  sufficiently  large,  therefore  L  —  •  y  >  C  •  y 

bounded  away  from  zero  and  the  minimum  above  can  be  computed  to 


min  e^®  —  we**®  —  ue“**®  =  Jl  —  |wol^  -  luol" 


=  \'-2(l-C“('-^))((dF+(pi)Zf) 


>  Ce. 


Combining  the  above  estimates  we  obtain  that 

eo>Ce-CM{kh)\ 

For  sufficiently  large  k  the  term  Cm  {khf  becomes  negligible  compared  to  e  = 
O  {k  {khy°),  resulting  in 

Co  ^  Ce  —  C  k  —  k 

completing  the  proof. 
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2.4.  A  stabilized  finite  element  method 

In  view  of  Theorem  2.6  it  is  clear  that  a  GFEM  with  no  pollution  term  in  the 
error  estimates  must  satisfy 

k  —  k  =  0. 

This  is  equivalent  to 

—  =  —2  cos  (kh) .  (2-14) 

N 

Of  course,  the  analysis  of  the  model  problem  (2.4)  guarantees  that  the  arising  so- 
called  “stabilized  finite  element  method”  (SFEM)  has  no  pollution  only  for  this 
special  example.  However,  two  of  the  authors  considered  the  more  general  problem 
(2.1)  and  showed  that  if  (2.14)  is  fulfilled,  it  is  always  possible  to  discretize  the 
boundary  conditions  and  the  inhomogeneous  right-hand  side  in  such  a  way  that 
the  GFEM  has  no  pollution.  The  details  are  in  the  following 


Theorem  2.7  (stabilized  finite  element  method).  Let  us  consider  the  prob¬ 
lem  (2.1).  Let  the  interval  (0, 1)  be  partitioned  using  the  grid  points  0  =  Xo  < 
xi  <  . ..  <  x„  =  1.  Here,  we  do  not  require  a  uniform  grid.  Let  the  nxn  system 
matrix  be  given  by 


_ 8in(fc(a;i4.i-Zj-i)) 

8in(k(*j+i-Xj))  8in(*(a:,— Zj_i)) 


1 

8in(fc|z<-z,|) 


0 


if  i  =j  <  n, 
if  i=j  =  n, 

if  \j  -  i\  =  1. 

otherwise 


and  the  mapping  by 


.  tan  (k^”‘  f(x)dx 

(QStob\  _  "  \  2  /  JXm-1  J  V  f - 

'  '•  2  tan  Xfn  3^171—1  (^m  ^m—l) 


The  corresponding  GFE-solution  u  is  defined  by 


Du  =  f. 


17 


while  u  denotes  the  exact  solution  of  (2.1). 

Under  the  assumption  that  the  right-hand  side  of  (2.1)  f  G'N}  (0, 1) ,  the  error 
estimate 

is  satisSed.  Obviously,  this  so-called  stabilized  Suite  element  method  (SFEM)  has 
no  pollution. 

Proof.  The  proof  of  this  Theorem  can  be  found  in  [2,  Theorem  2.3]. 

Remark  2.  In  the  case  of  a  uniform  mesh,  for  the  SFEM,  the  quotient  D2/N 
satisfies 

— 2  cos  (kh) . 


3.  The  GFEM  to  the  Helmholtz  equation  in  two  dimen¬ 
sions 

In  this  section  we  consider  the  Helmholtz  equation  on  a  square  of  side  length  2L 

-  Alt  -  k^u  —  f  in  =  {—L,  Vf  (3.1) 


with  boundary  conditions 


asiku  +  =  g»  on  T*  for  s  G  {1, 2, 3, 4} .  (3.2) 

on 

The  boundary  pieces  F,  are  depicted  in  figure  (3.1)  where  the  symbol  ^  means 
“outside  normal  derivative” . 

Let  fl  be  partitioned  into  squares  of  side  length  h  =  forming  the  finite 
element  grid  r.  The  space  of  finite  element  functions  5/,  is  defined  by  continuous, 
piecewise-bilinear  functions  corresponding  to  the  grid  r.  The  grid  points  are 
defined  by  Xs,t  ■=  h{t-l,s-  1)^.  The  local  basis  functions  are  denoted  by 
and  satisfy  the  relation 


^8,t  (^S'jt') 


1  if  (s,  t)  =  (s',  If) 
0  otherwise. 
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Figure  3.1:  Domain  Cl  with  boundary  and  normal  directions 

A  finite  element  function  u  £  Sh,  can  be  identified  with  a  vector 
the  basis  representation 

« (^)  =  H  ^s,t(l>s,t  (a^)  •  (3-3) 

s,t=l 

The  GFEM  is  a  method  which  defines  an  x  matrix  A  and  a  linear  mapping 
Q  which  maps  the  right-hand  sides  of  (3. 1,3.2)  onto  the  vector  of  the  right-hand 
side  b,  i.e. 

b  :=  Q  {g,  f)  ■ 

The  solution  of  the  linear  system 

Au  =  b 

is  identified  with  a  finite  element  function  by  (3.3)  which  serves  as  an  approxima¬ 
tion  of  the  exact  solution  of  (3.1, 3.2). 

The  dimension  of  the  system  matrix  is  thus,  each  grid  point  Xa,t  can  be 
associated  with  one  equation.  If  the  grid  point  Xa,t  is  an  interior  grid  point,  i.e. 
not  lying  on  the  boundary,  then  the  corresponding  equation  can  be  written  in  the 
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lorm  ,  j 

U«_i,t+i  +Ag\tUg^t+l  Ua+i,t+i 

+  +^“;tUg,t  +>lJ|°Ua+i,t 

+  ’  ^Ua_i,t_i  +-<42, ’t  ^Ug,t_i  +^i’t  ^Ug+i,t_i  =  6a, t- 

If  the  grid  point  is  lying  on  the  boundary  F,  then  those  elements  of  .<4^^  which 
are  multiplied  with  values  of  Up^q  with  h{q—  l,p  —  1)^  ^  have  to  be  set  to  zero. 
This  means  that  the  system  matrix  of  a  GFEM  is  set  up  by  defining  all  interior 
“stencils”  , ,  „ ,  , ,  _ 


A  interior 
■^s,t 


^7}’' 


j-1,0  40,0  41,0 

• -  ■'^*.4  -^Sft  -^Syt 

A-h-'^  Al>7^ 

•  *'Bft  Syt  OfZ 


all  edge  stencils 


j^edge 

Sytr 


A7l’^ 


and  all  corner  stencils 


A  corner  _ 

^8yt  ^ 


A7,l'" 

A7l'° 

A7}'-^ 


A7r 

A7}’° 


The  quality  of  our  GFEM  approximation  depends  on  the  coefficients  A[’t  and  the 
definition  of  the  right-hand  side  vector.  We  impose  the  following  restrictions  on 
the  stencils  -<4a,t  and  mapping  Q. 

A1  The  interior  stencils  which  depend  on  k  and  h,  satisfy  the  following 

symmetry  condition 

A2  Ai  A2 

^interior  _  (3.4) 

A2  Ai  A2 


and  have  the  same  values  for  all  s  and  t. 
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A2  Let  us  assume  that  the  right-hand  side  /  of  (3.1)  is  zero.  Then,  the  lin¬ 
ear  operator  Q  is  local  in  such  a  way  that  if  Xg^t  Is  a  grid  point  satisfying 
dist  {xg,t,  r)  >  2h  then  the  corresponding  entry  of  the  right-hand  side 
vector  is  zero. 

A3  We  assume  that  the  interior  stencils  of  the  finite  element  matrices 
can  be  expanded  as  a  Taylor  series  of  the  form 

(i)  = 

(ii) 

(iii)  = 

with  a  =  kh  and  independent  of  k  and  h  for  alH  G  {0, 1, 2} ,  m  G 

{0,1,2,...}. 


A4  We  assume  that  the  principal  part  of  A,  i.e. 


Aprincipal  • — 


(^2)0  (Al)o 
(Al)o  (Ao)o 
(^2)0  (Al)o 


(^2)0 

(■^1)0 

UzSo  . 


is  an  approximation  of  the  principal  part  cq  (u,  v)  =  (Vu,  Vv)dx  oi  con¬ 
sistency  order  2,  implying 


(Ao)o  >  0, 


(■^0)0  +  4  ((Ai)o  +  (^2)0)  —  0,  (3.5) 


-(A0o-2(A2)o  =  1. 

These  restrictions  are  very  natural  considering  linear  finite  elements.  Some 
comments  are  given  in  the  following 

Remark  3.  Condition  A1  reflects  to  the  rotational  and  translational  symmetry 
of  the  Helmholtz  equation  and  the  mesh  r. 

Condition  A2  reSects  the  fact  that  the  discretization  of  the  boundary  condi¬ 
tions  has  local  influence  to  the  right-band  side  vector. 

Condition  A3  reflects  the  fact  that  the  Laplacian  and  the  identity  operator  are 
of  even  order. 

Condition  A4  is  the  usual  consistency  condition  if  "—A”  is  discretized  by  linear 
elements. 
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4.  Approximation  of  the  Helmholtz  equation  by  the  GFEM 


In  this  section  we  will  investigate  the  dependency  of  the  accuracy  of  the  GFEM 
on  the  matrix  stencils.  To  measure  the  accuracy,  or  the  difference  between  the 
GFE  solution  and  the  exact  solution,  we  introduce  a  weighted  £^-norm  defined 

ll.'lli  :=  / 

Jn  1+  ||a:|| 

This  weighted  norm  refiects  the  fact  that  in  this  paper  our  aim  is  not  the  modeling 
of  the  DtN  boundary  conditions  and  their  discretization  but  to  discretize  the 
Helmholtz  operator  in  the  interior  of  the  domain  in  an  optimal  way. 

To  specify  the  quality  of  our  GFE  discretization  we  will  use  some  tools  of  the 
theory  of  the  (integral)  Fourier  transform  and  the  discrete  Fourier  transform.  We 
give  here  only  a  short  summary  of  the  theory  presented  in  [2].  The  symbol  of  the 
Helmholtz  operator  is  given  by 

where  ^  and  1|^||^  :=  +  ^1- 

In  Condition  A1  of  the  previous  section  we  assumed  that  the  interior  stencils 
of  the  GFE-matrix  have  constant  coefficients,  i.e. 


j^interior  _ 

Sft 


A2  Ai  A2 
Ai  Aq  Ai 
A2  Ai  A2 


where  At  only  depends  on  k  and  h.  The  discrete  symbol  of  the  corresponding 
difference  operator  is  given  by 

(T stencil  (0  ==  +  2Ai  (cOS  6  +  COS  6)  +  4^2  COS  ^  COS  6- 

Let  Afkh  be  defined  as  the  scaled  roots  of  the  operator  symbol  a: 

:=  {«  €  I  a  (ft-'e)  =  o} ; 

in  other  words  Afkh  is  a  circle  centered  at  the  origin  with  radius  kh. 

Further,  let  NstencU  be  defined  as  those  roots  of  (JatencU  lyiag  in  the  square 
-  e,  +  e)  x  {-kh  -t,kh->r  e)  where  e  >  0  has  to  be  chosen  sufficiently 
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small  such  that  Mkh  is  a  simple  connected  line.  The  maximal  distance  between 
these  curves  defined  by 


V  {stencil)  :=  V  := 


(4.1) 


can  be  considered  as  a  measure  of  the  approximation  quality  of  the  GFEM  for  the 
Helmholtz  equation  in  the  interior  by  the  GFEM  difference  operator.  The  details 
can  be  found  in  the  following 


Theorem  4.1. 
problem 


Let  be  the  interior  stencil  of  a  GFEM  to  solve  the  Helmholtz 

—  Alt  —  k^u  =  f  in  Ql  '•=  {—L,  L)  x  {—L,  L)  (4.2) 


with  boundary  conditions 

Bu^g  (4.3) 

which  should  imply  existence  and  uniqueness  of  the  solution. 

a)  Then  there  exists  L  <  oo,  a  right-hand  side  /,  and  boundary  data  g  such 
that  the  error  of  the  GFE-solution  u/e  can  be  estimated  from  below  by 

On  the  other  hand  there  exists  a  function  Uopt  in  the  Snite  element  space  which 
satisfies 

||u  -  UoptW-  <  C2k^h'^, 

where  the  constants  Ci  and  Cz  are  independent  of  k  and  h. 

b)  The  function  V  {stendl)  can  be  expanded  accordingly: 

V  {stendl)  =  rjj,  {kh)^^'^^  +  O  {kh)^°'^^ 

with  constants  ri  depending  on  the  stencil  coefficients  {At)^  for  t  e  {0,1,2}  ,m  £ 
No  (cf  Condition  A3)  but  not  on  k  and  h.  Further,  1<  lo  <  oo  and  the  coefficient 

rio  7^  0-  V  M  II  I. 

c)  Asymptotically  (i.e.  for  sufficiently  small  kh)  the  error  yu  -  ii/el|_  can  be 
estimated  by 

||lt  —  U/e||_  >  CiCftencil'^  (^^)  **  • 

Consequently,  if  kh  is  small  enough,  the  error  of  the  best  approximation  could  be 
small,  but  for  sufficiently  large  k  the  error  of  the  GFE-solution  could  be  arbitrarily 

large. 
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Proof.  The  proof  of  this  Theorem  is  given  in  [2,  chap.  3].  ■ 

The  theorem  can  be  interpreted  as  follows.  Let  a  GFEM  be  given.  Then  the 
following  situation  can  arise  for  given  right-hand  side  /  and  boundary  data  g. 
We  want  to  compute  the  solution  of  (4.2, 4.3)  with  an  accuracy  of  e.  Then  by  the 
approximation  property  of  our  finite  element  space  Sfi  we  know  that  if  the  number 
of  elements  is  larger  than  no  there  exists  a  function  Uopt  £  <5/,  satisfying 

II M  -  UoptW  <  e.  (4.4) 

Then,  the  number  of  elements  ugfem  to  guarantee  that  the  GFE-solution  also 
fulfills  (4.4)  has  the  property  that  the  ratio  ncFEMlna  behaves  asymptotically 
like  „  , 

no 

Obviously  the  ratio  goes  to  infinity  with  increasing  k.  We  conclude  that 

the  function  V  (stencil)  is  a  well-suited  measure  of  the  approximation  quality  of 
the  GFE  discretization  to  solutions  for  the  Helmholtz  equation. 

5.  A  generalized  finite  element  method  having  minimal 
pollution 

Based  on  the  theoretical  results  of  the  previous  section  we  are  now  able  to  con¬ 
struct  a  GFEM  having  minimal  pollution.  To  be  more  concrete  we  will  define  an 
interior  stencil  such  that  V  (stencil)  is  asymptotically  minimal. 

By  Assumption  A3  of  Section  2,  Aq  0,  therefore  the  quotients  Oi  and  02  are 
well-defined  by 

“2=4^. 

Using  A4,  Cl  and  02  can  be  expanded  as 

00 

m=0 

and 

02 = i:  (kh)^  ■  (5.1) 

m=0 
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Note  that  condition  (3.5)  implies  that 


(5.2) 


1  +  Ao  +  ifl  ~  0 
Ao  +  2io  ^  0- 

In  view  of  the  definition  of  X>  (stencil)  (cf.  4.1)  we  introduce  the  function  dist  (-tt ,  tt) 

Ro  by 


/  cosp  \ 

V  sin/3  ) 

The  function  dist  can  be  expressed  explicitly  by  the  expansion 


dist  (^)  :=  min  fch 


dist (P)  := 


E  rm  (0)  (kh) 


2m+l 


m=l 


Before  we  define  the  coefficients  =  Vm  (0)  we  introduce  Kn,  t„  and  P2n,m  by 


f— coe' 

Kn  = 


(2n)! 


p2n,m  —  ^ 


Tfx  =  ^2n)!  *  "^)n  ^m=0  ^mTn—ri 

4,m  if  n  =  0, 


r  -krif . .  .-kr 


otherwise 


\2n-fold  convolutum/  . 


with  dn,m  denoting  the  Kronecker  delta. 

Formally  we  set  Tq  =  1.  Then,  all  other  coefficients  (P)  are  given  by  the 
condition 


n=0  m=0  \  Z  / 

We  state  that  condition  (5.3)  can  be  written  in  the  form 

The  abbreviation  l.o.t.  denotes  the  remaining  sum  of  (5.3)  containing  only  func¬ 
tions  rj  with  j  <  /  - 1.  In  view  of  (5.2),  relation  (5.4)  serves  as  a  recursive  formula 
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for  the  functions  rj.  The  proofs  and  development  of  these  formulae  can  be  found 
in  [2,  Appendix]. 

In  the  mentioned  paper  it  was  further  proved  that,  for  bounded  kh  <  ao) 


X>  {stendl) 


max  dist  (P)  <  C 


max 

—7r<P<7r 


a(/5)  > 


(5.5) 


where  the  largest  possible  value  of  Iq  is  3. 

To  summarize  this  section,  we  state  that  the  quality  measure  "D  [stencil)  can 
be  explicitly  computed  by  formulae  (5.3)  and  (5.5)  for  each  generalized  finite 
element  method. 

By  the  consideration  above  it  is  clear  that  an  asymptotically  optimal  interior 
stencil  has  to  be  designed  such  that  ri  [0)  =  [0)  =  0.  According  to  condition 

A3  (cf.  Section  3),  we  had  assumed  that  the  interior  stencil  Ainterior  can  be  written 
in  the  form 


A2  Ai  A2 

00  00 

Am,2  Ato,1  Afnfl 

j^interior  _ 

Ai  Aq  Ai 

=  •£  =  £  (kh)^ 

Afn^l  Affifi  Affi,l 

A2  Ai  A2 

m=0  m=0 

Am, 2  ^m,l  ^m,2  . 

(5.6) 

where  the  interior  stencils  in  the  expansion  above  are  independent  of  k 

and  /i,i.e.,  Am,t  are  in  general  complex  numbers.  In  the  following  we  will  use 
a  :=  kh. 

We  define  the  interior  stencil  A^pf%f  of  the  quasi- stabilized  FEM  (QSFEM) 
by 

Ao  =  4 


A  _  o _ ci(a)8i(o)-C2(a)s2(a)  ,f. 

^C2(a)a2(a)(ci(a)+«i(a))-ci(Q)»i(a)(c2(a)+«2(a))  v  ‘  I 


Ao  — 


C2(a)+»2(a)-ci(o)-gi(a) 


C2(a)s2(a)(ci(a)+«i(a))-ci(a)«i(o)(c2(a)+»2(a))  ’ 

while  the  auxiliary  functions  Ci,Si,C2,  and  S2  are  defined  by 
Cl  (a)  :=  cos  (a  cos  Si  (a)  :=  cos  (a  sin 

C2  (a)  :=  cos  (a  cos  f|)  S2  (a)  :=  cos  (a  sin  f|)  . 

The  function  Ai  and  A2  be  expanded  according  to  (5.6)  because  cos  is  an  even 
function.  The  first  terms  of  this  expansion,  i.e.  the  stencils  for  m  <  4, 
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are  given  by 


I  interior  _ 


1  _4 

1  ® 

--  4 

5  _  4 

"5  5 


I  i  _| 


i  interior  _  _ 


Ainterior  _  _ 

wHo  — 


76313  473849  76313 

'^§W  473849  _:!W 

■  22500000  45000000  22500000 


j^interior  _  _ 


826713271 


5094901033 

’2376000000000 

0 

5094901033 


826713271 


mmr 


-4  -  23^^g^^§^^j00  .^901033 _ , 

L  1188000000000  2376000000000  1188000000000  J 

The  properties  of  this  stencil  can  be  found  in  the  following 

Theorem  5.1.  Let  deSned  by  (5.7).  This  stencil  has  the  property 

that  the  functions  Tm  (P)  deSned  by  (5.3, 5.4)  satisfy 

ri  iP)  =  r2  iP)  =  0 

/  m  C06(8^2 
^3  KP)  774144  • 

The  quantity  V  can  be  estimated  by 

D  (4"^:™)  <  1-3  ■  10-*  (*A)’  +  o  {{khf)  • 

Proof.  The  proof  is  purely  technical  but  simple,  therefore  we  give  here  only 
an  outline  of  it. 

First  one  has  to  expand  the  quotients  4  and  4^^^^  according 

V  QSFEM rtctPDKM  ^ 


(  Ainterior  \ 
\^QSFEM)q 


to  (5.1)  to  obtain  the  coefficients  Xm  and  The  proof  then  is  completed  by 
computing  the  functions  (P)  using  the  recursive  formula  (5.4).  ■ 
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6.  Numerical  Results 

In  this  section  we  will  present  the  results  of  an  implementation  of  different  GFEM 
to  approximate  the  following  problem 

-Au-k^u  =  0mQ:=  (0, 1)  x  (0, 1)  (6.1) 

with  boundary  conditions 

iku  +  —  —  g  onT  :=  dCl.  (6.2) 

on 

The  function  g  depends  on  the  parameter  6  and  is  given  by 

'  i{k-k2)e^^^^  ifxGri:=(0,l)x(0,0), 

i  {k  +  ki)  e‘(*‘+****)  if  a;  G  Ta  :=  (1, 1)  x  (0, 1) , 

g  {x)  :=  ^ 

i  {k  +  ki)  if  X  G  Ts  :=  (0, 1)  x  (1, 1) , 

i{k-  ki)  if  X  G  r4  :=  (0, 0)  x  (0, 1) 

with  (fci  yki)  =  k  (cos  6,  sin  6) . 

The  exact  solution  of  this  problem  is 

Hex  (a:)  :=  +****>. 

6.1.  Discretization  techniques  for  the  Helmholtz  equation 

We  discretize  the  domain  by  squares  of  side  length  h  =  and  use  bilinear 
elements.  Consequently  the  system  matrix  has  dimension  n^.  We  implemented 
the  following  three  discretization  methods. 

1  Galerkin  Finite  Element  Method 

Writing  (6.1)  in  a  variational  formulation  and  incorporating  the  boundary 
conditions  results  in  the  following  problem. 

Find  u€'H^  (Q.)  such  that 

J  VuVu  —  k^uvdx  +  uvdTx  =  gvdTx 
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is  satisfied  for  all  v  G  (f^).  Replacing  the  infinite-dimensional  space  {Q)  by 
the  finite  element  space  Sh  of  bilinear  elements  and  introducing  the  usual  local 
basis  results  in  a  system  of  linear  equations  for  the  coefficients  u  of  the  basis 
representation  (3.3).  This  method  can  also  be  described  in  terms  of  stencils.  The 
interior  stencil  for  the  Galerkin  method  is  given  by 


Ainterior 

■^gal 


■■=  t 

.  3  3  3  - 


i  ?  i 

36  9  36 


A  similar  computation  as  explained  in  the  proof  of  Theorem  5.1  3delds  the  fol¬ 
lowing  estimate  for  the  quantity  V  {stendl)  which  describes  the  approximation 
quality  of  this  method  for  the  Helmholtz  equation. 

®  KJ-'-)  =  (khf  +  O  {(khf)  =  ^  +  O  Hkhf)  . 

2  Generalized  Least  Squares  Finite  Element  method  (GLS-FEM) 

In  [9],  Thompson  and  Pinsky  have  generalized  the  GLS-FEM,  originally  intro¬ 
duced  by  Harari  and  Hughes  in  [4]  for  a  one-dimensional  model  problem,  to  the 
two  space  dimensions.  Applied  to  problem  (6.1)  discretized  by  bilinear  elements 
this  method  can  be  written  in  the  form: 

Find  Sh  such  that 

j  VuVu  -  Tuvdx  +  j  uvdTx  =  gvdTx 

is  satisfied  for  all  u  G  5/,.  The  parameter  r  =  r  (fc,  h)  is  given  by 

4  -  cos  (kh  cos  f )  -  cos  (kh  sin  f)  -  2  cos  (kh  cos  f )  cos  (kh  sin  f ) 

^  ^  ^2  +  cos  (kh cos  l])  (2  -I-  cos  (kh sin  h^ 

The  interior  stencil  for  the  GLS-FEM  can  be  written  in  the  form 

r  _i  _i  _i  T  r  -L  1  -L  1 

333  „  36  9  36 

Ainterior  _!  8  —h^rik  h^  ill 

■^GLS-FEM  ■—  5  ’  f  ?  ? 

.  “3  “3  ”3  J  L  ^  9  ^  . 

To  compare  this  method  with  the  standard  Galerkin-FEM  we  state  that 


-h\  =  -  {khf  -  j-  (hk)^  +  O  ({khf) 
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which  means  that  the  GLS-FEM  for  bilinear  elements  is  a  higher  order  modifi¬ 
cation  of  the  interior  stencil  of  the  Galerkin  FEM.  The  approximation  quality  of 
this  method  to  the  Helmholtz  equation  in  terms  of  V  {stendl)  is  given  by 

_m^^  +  O  {{khf)  =  +  O  ((kh)^) 

and  hence  the  pollution  is  essentially  the  same  as  for  the  Galerkin  FEM. 

3.  Quasi-Stabilized  Finite  Element  Method  (QSFEM) 

The  interior  stencil  of  this  method  was  already  presented  in  Section  5.  Here 
we  describe  only  the  modeling  of  the  boundary  conditions  and  the  assembling  of 
the  vector  of  the  right-hand  side.  This  is  done  analogously  to  the  finite  difference 
method  by  replacing  the  normal  derivatives  by  a  difference  formula  centered  at  the 
edge  points  and  then  eliminating  the  fictitious  points.  This  technique  is  described 
together  with  an  error  analysis  in  [3,  Section  4.7.2]. 

We  recall  that  the  approximation  quality  of  this  method  was  proved  to  be 

p  Ki?r«)  =  (*'■)'.:“<?<,  +0 

6.2.  Numerical  evaluation  of  the  GFE  methods  for  the  Helmholtz  equa¬ 
tion 

In  order  to  measure  the  accuracy  and  compare  the  presented  methods  we  have  to 
introduce  suitable  norms.  A  natural  measure  of  the  approximation  error  of  (6.1) 
is  the  energy  norm,  i.e.  the  usual  ?{^-seminorm 

_  li^ez  ^/elli 

with 

ll^ll?  :=J^Vu{x)-Vu  (x)  dx. 

Alternatively  we  will  measure  the  accuracy  of  the  solution  in  the  £^-norm 

_  ll^ex  ~  ^/ello 

with  . 

Iloilo  /  u{x)u{x)dx. 

J 
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Recalling  the  one-dimensional  results  (cf.  Section  2)  we  expect  that  the  error 
of  the  best  approximation  behaves  like 


Ugx  I 


<  Chk, 


and  - 


In  the  one-dimensional  case,  the  error  of  the  Galerkin-FEM  could  be  estimated 

ei<C  (kh)  -I-  Ck^h  {kh) 


eQ<C{hk)^  +  Ck{khf.  (6.4) 

Obviously,  the  pollution  of  the  ?^^-error  becomes  negligible  if  k^h  is  small,  but  in 
the  pre-asymptotic  range  we  expect  that  the  Galerkin  solution  differs  substantially 
from  the  best  approximation. 

For  the  £^-error  the  pollution  term  is  the  dominant  term  in  the  error  estimate 
for  all  values  of  h.  We  expect  that,  with  increasing  value  of  k,  the  distance  of  the 
graph  of  the  Galerkin  error  from  the  best  approximation  error  increases. 

The  error  of  the  GFE  solution  depends  on  the  direction  6  of  the  wave  vector 
k,  where 

k  =  {ki,k2]^  =  k  (cos  6,  sin  6) . 

It  turns  out  that  for  special  values  of  9  the  GFE  solutions  effectively  coincide 
with  the  best  approximation,  where  for  other  values  of  9  the  GFE  solution  differs 
substantially  from  the  best  approximation.  Unfortunately,  the  direction  9  is  not 
known  a  priori.  In  order  to  guarantee  that  the  error  of  the  GFE  solution  is  under 
control  (i.e.  the  solution  is  reliable)  one  has  to  assume  that  the  error  is  sufficiently 
small  even  if  the  solution  would  correspond  to  the  value  of  9  with  the  largest  error. 
Therefore  for  j  G  {0, 1},  we  have  computed  the  quantity 

=  max  e,-  (a;) . 

Here,  ej  (cj)  denotes  the  error  of  the  GFE  solution  of  problem  (6.1, 6.2),  where  the 
parameter  9  for  the  right-hand  side  g  is  chosen  as  u>.  The  exact  solution  in  this 
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case  is  given  by 

We  approximate  the  maximum  above  by  choosing  the  set  of  a;-values 

_  TT  ttSthttStt') 


^  f-  TT  TT  StT  TT  SttI 

t^’Te’ 8’16’4’Tj 


Due  to  the  rotational  symmetry  of  our  problem  we  know  that  the  error  corre¬ 
sponding  to  a  direction  6  is  the  same  as  for  0  +  for  integers  m.  The  maximum 
above  is  approximated  by 


We  restrict  the  step  size  h  to 


«  maxe,-  (a;) . 

^  oigH 


hk  < 

-  2 


This  assumption  is  natural  because  it  guarantees  at  least  four  elements  per  wave 
length  (see  Fig.  6.1). 


Approximation  of  sin(x)  by  linear  elements  witb  4  elements  per  wave  length 

1  I ’ 

0.8  -  /y'  \\  8in(x) - 

//  \\  linear  interpolant - . 


2  3  4  5  6  7 


Figure  6.1:  Approximation  of  sin(x)  by  the  piecewise  linear  interpolant  with  four 
elements  per  wave  length. 

For  larger  values  of  h  the  error  of  the  best  approximation  would  be  too  large 
for  practical  applications. 

Figure  6.2  and  6.5  show  the  H^-evior,  resp.  £2-error  of  the  three  discretization 
schemes  and  of  the  best  approximation  for  k  =  30, 100, 150.  The  plots  consists 
of  groups  of  four  lines,  where  the  lowest  line  always  corresponds  to  the  best  ap¬ 
proximation,  the  next  one  corresponds  to  the  QSFEM,  the  third  belongs  to  the 
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GLS-FEM  and  the  highest  line  always  corresponds  to  the  Galerkin  solution.  We 
see  that  the  error  of  the  Galerkin  solution  behaves  as  expected.  The  V}-etTOT 
differs  substantially  from  the  best  approximation  in  the  pre-asymptotic  range, 
while  this  range  increases  for  larger  values  of  k.  On  the  other  hand  the  pollution 
becomes  negligible  if  h  is  sufficiently  small,  since.  k‘^h  is  small.  The  situation  for 
the  £^-norm  is  different.  Here,  the  pollution  does  not  vanish  as  h  — >  0.  The  lines 
becomes  parallel  (in  our  log-log  plot)  for  small  values  of  h  whereas  the  distance 
from  the  best  approximation  increases  with  higher  wave  number  k  for  all  h  in 
accordance  with  the  theoretical  estimate  (6.4). 

The  improvement  of  the  QSFEM  is  obvious.  The  size  and  range  of  the  pol¬ 
lution  even  for  A:  =  150  is  very  small  and  remains  nearly  constant  for  different 
values  of  k.  For  the  £^-norm  this  behavior  can  be  observed  from  the  constant 
(w.r.t.  the  wave  number  k)  distance  between  the  graph  of  the  best  approximation 
and  the  QSFE-solution. 

The  GLS-FEM  shows  nearly  no  improvement  over  the  Galerkin  FEM  for  rel¬ 
atively  large  values  of  hk  ~  7r/2,  while  the  pollution  decreases  faster  with  respect 
to  h  than  for  the  Galerkin  FEM.  These  numerical  results  are  in  accordance  with 
the  different  sizes  of  the  theoretical  quality  measure  V  {stendl)  computed  above. 
The  pictures  (6.2-6.5)  can  also  be  interpreted  as  follows.  Let  us  assume  that  we 
want  to  approximate  the  solution  of  our  model  problem  with  an  relative  error  of 
€  in  the  £^-norm.  Then  we  can  ask  how  many  degrees  of  freedom  (DOF)  are 
necessary  to  get  this  accuracy.  The  following  table  shows  this  dependency  for 
some  values  of  e. 


DOF  necessary  to  obtain  an  accuracy  of  e  in  the  £‘^-norm 

=  100 

k  =  m 

e 
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ESSHI 

TOiia 

As  mentioned  above  the  error  of  the  GFE-solution  depends  significantly  on  the 
wave  direction  9.  In  the  following  plots  we  illustrate  the  pollution  effect  depen¬ 
dent  of  the  parameter  6.  For  constant  kh  we  know  that  the  error  of  the  best 
approximation  in  the  T^Lnorm  is  of  order  kh.  We  have  chosen  kh  =  1.5, 0.7, 0.3, 
where  kh  =  1.5  «  7r/2  corresponds  to  four  elements  per  wave  length,  kh  =  0.7 
to  eight  and  kh  =  0.3  to  16  elements  per  wave  length.  In  Figures  6.6-6.14,  we 
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have  plotted  the  difference  of  the  GFE-errors  from  the  best  approximation  error. 
We  see  that  the  Galerkin  FEM  is  relatively  independent  of  the  parameter  B  but 
the  difference  from  the  best  approximation  is  significantly  for  all  considered  cases. 
The  difference  of  the  GLS-FE  solution  from  the  best  approximation  is  practically 
zero  for  0  =  I ,  where  for  0  =  0,  f ,  f  it  is  especially  for  kh  =  7r/2  practically  as 
bad  as  the  Galerkin  FEM.  The  scaling  of  the  axes  of  the  GLS-FEM  plots  and  the 
QSFEM  plots  is  always  the  same.  We  see  that  for  the  considered  range  of  k  the 
QSFEM  has  practically  no  pollution.  The  difference  from  the  best  approximation 
is  small  for  all  considered  values  of  6  and  is  fairly  steady  with  the  wave  number 
k.  In  contrast  to,  e.g.  6  =  0,  the  difference  of  the  GLS-FE-solution  from  the  best 
approximation  increases  with  higher  wave  number  k  which  means  for  this  method 
it  is  not  sufficient  to  restrict  the  quantity  kh  to  get  a  small  error.  Figures  6.15- 
6.23  shows  the  dependency  of  the  pollution  on  kh  and  6  in  the  £^-norm.  One  can 
observe  that  the  pollution  of  the  QSFEM  is  even  smaller  than  in  the  ^^-norm, 
where  the  behavior  of  the  GLS-FEM  and  Galerkin  FEM  is  quite  similar. 

6.3.  Conclusions  and  Remarks 

In  this  subsection  we  summarize  and  comment  on  the  numerical  results  and  relate 
them  to  the  theory  presented  in  the  previous  sections.  In  Section  2  we  had  ana¬ 
lyzed  a  one-dimensional  example  which  shows  the  typical  pollution  effect  for  the 
GFE-approximations.  We  have  seen  that  if  we  are  able  to  design  a  GFEM  which 
has  no  pollution  for  this  example  it  will  also  have  no  pollution  for  more  general 
cases  (cf.  Theorem  2.7)  if  we  treat  the  boundary  conditions  and  an  inhomogenous 
right-hand  side  in  an  appropriate  way.  On  the  other  hand  we  have  proven  that  if 
the  discrete  wave  number  and  the  exact  one  are  not  the  same,  then  the  pollution 
cannot  be  avoided  by  any  treatment  of  the  boundary  conditions. 

In  the  two-dimensioned  case  we  have  seen  that  a  necessary  requirement  to  get  a 
small  pollution  is  to  minimize  the  maximal  distance  V  of  the  zeros  of  the  discrete 
symbol  and  of  the  operator  symbol.  We  were  able  to  design  a  GFEM  called 
QSFEM  such  that  V  is  minimal.  On  the  other  hand  we  have  not  investigated  a 
special  treatment  of  the  boundary  conditions.  Nevertheless,  from  the  numerical 
results  presented  above  it  follows  that  even  a  straightforward  treatment  of  the 
boundary  condition  (centered  differences)  results  in  a  GFEM  with  a  negligible 
pollution  effect  for  moderate  wave  numbers.  We  conclude  that  on  the  one  hand  the 
pollution  effect  for  the  FE  treatment  of  the  Helmholtz  equation  is  not  avoidable 
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in  principle,  while  on  the  other  hand  it  is  possible  to  design  the  QSFEM  having 
a  pollution  which  is  not  visible  in  practice. 
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relative  error  In  the  H1  -seminorm  relative  enor  In  the  H1  -seminorm 


Figure  6.2:  Relative  error  in  the  T/^-seminorm  for  k  =  30, 100 


Figure  6.3:  Relative  error  in  i7^-seminorm  for  k  =  150 
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relative  error  in  the  L2-norm  relative  error  In  the  L2-nortn 


Figure  6.4:  Relative  error  in  the  L^-norm  for  k  =  30, 100 


Figure  6.5:  Relative  error  in  the  L^-norm  for  k  =  150 
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el  (GLS)-el  (best  approximation)  ^  el(Galerkln-FEM)-el  (best  approximation) 


6.6:  Dependency  of  the  H^-exTor  on  the  angle  6  for  kh  =  1.5  for  the 

Galerkin-FEM 


Figure  6.7:  Dependency  of  the  H^-exrox  on  the  angle  6  for  kh  —  1.5  for  the 

GLS-FEM 
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Figure  6.8:  Dependency  of  the  H^-enoT  on  the  angle  9  for  kh  =  1.5  for  the 

QSFEM 


Figure  6.9:  Dependency  of  the  H^-erroT  on  the  angle  6  for  kh  —  0.7  for  the 

Galerkin-FEM 
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Figure  6.10:  Dependency  of  the  H^-evvor  on  the  angle  6  for  kh  —  0.7  for  the 

GLS-FEM 
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Figure  6.11:  Dependency  of  the  i/^-error  on  the  angle  6  for  kh  =  0.7  for  the 

QSFEM 
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el  (GLS)-e1  (best  approximation)  ^  o1(Galerkln-FEM)-o1  (best  approximation) 


approximation)  ^  e1(QSFEM)-e1  (bast  approximation) 
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e  6.14:  Dependency  of  the  H^-erroT  on  the  angle  6  for  kh  =  0.3  for  the 

QSFEM 
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Figure  6.15:  Dependency  of  the  L^-error  on  the  angle  9  for  kh  =  1.5  for  the 

Galerkin-FEM 


42 


Figure  6.16:  Dependency  of  the  L^-error  on  the  angle  9  for  kh  =  1.5  for  the 

GLS-FEM 
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Figure  6.17:  Dependency  of  the  L^-error  on  the  angle  9  for  kh  =  1.5  for  the 

QSFEM 
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Figure  6.18:  Dependency  of  the  L^-error  on  the  angle  6  for  kh  =  0.7  for  the 

Galerkin-FEM 


Figure  6.19:  Dependency  of  the  L^-error  on  the  angle  6  for  kh  =  0.7  for  the 

GLS-FEM 
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eO(GaIerkin-FEM)-eO(be8t  approximation)  C  eO(QSFEM)-aO(best  approximation) 


■e  6.20:  Dependency  of  the  L^-error  on  the  angle  9  for  kh  —  0.7  for  the 

QSFEM 


Figure  6.21:  Dependency  of  the  L^-error  on  the  angle  6  for  kh  =  0.3  for  the 

Galerkin-FEM 
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eO(QSFEM)-eO(beat  approximation) 


wave  direction  theta 


Figure  6.22:  Dependency  of  the  L^-error  on  the  angle  6  for  kh  —  0.3  for  the 

GLS-FEM 


Figure  6.23:  Dependency  of  the  L^-error  on  the  angle  6  for  kh  =  0.3  for  the 

QSFEM 
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